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I. INTRODUCTION 

To understand hadronic and nuclear scattering at the 
high collision energies reached at the LHC, one needs a 
good picture of high energy behavior of QCD. Much of 
the phenomenological work in the field is done within 
the Color Glass Condensate (CGC) framework (for re- 
views see e.g. [1, 2]). In this picture, the most convenient 
choice for the basic degrees of freedom of a hadronic state 
are Wilson lines (path ordered exponentials of the color 
field), which describe the eikonal propagation of a high 
energy probe. In the CGC, these are stochastic vari- 
ables drawn from an energy-dependent probability dis- 
tribution. This energy dependence is described by the 
JIMWLK renormalization group equation (see Ref. [3- 
5]), which can be derived in weak coupling QCD by suc- 
cessively integrating additional gluonic degrees of free- 
dom into the Wilson lines. 

With a solution of the JIMWLK equation one can cal- 
culate expectation values for various operators formed 
from the Wilson lines, describing different scattering pro- 
cesses. For example the simplest one, the dipole, (see 
Eq. (10)) is related to the inclusive DIS cross section and 
to the unintegrated gluon distribution, needed for sin- 
gle inclusive particle multiplicities. In phenomenological 
work, the dipole can most conveniently be obtained from 
the Balitsky-Kovchegov (BK) equation [6-8] which, al- 
though the original derivation predates JIMWLK, can 
now be viewed as its mean field approximation. Multi- 
particle correlation observables, however, require an un- 
derstanding of the JIMWLK evolution beyond the BK 
approximation. For example, for a calculation of the az- 
imuthal angle dependence of two-particle correlations [9- 
13] one needs the expectation values of four and six point 
functions such as the quadrupole, Eq. (28). Another 
example are long range rapidity correlations, e.g. the 
"ridge" , which is related to long range correlations in the 
color fields of the colliding hadrons [14-18]. 

The inclusion of the running coupling constant in 
the BK equation (rcBK) [19-21] has been an essential 
theoretical improvement for phenomenological applica- 
tions [22-26]. The running of the QCD coupling slows 



down the evolution (with x or equivalcntly In yfs) signif- 
icantly from the fixed copupling case, bringing the speed 
to a rough agreement with HERA measurements. Since 
the full NLO version of the equation [27] has not been ex- 
tensively applied to phenomenology (see however recent 
work in Ref. [28]), rcBK is currently the state of the art 
in phenomenological applications. 

The running coupling resums a certain subset of all 
the NLO corrections to the evolution equations. Which 
subset precisely one chooses to resum is not uniquely de- 
fined, but different choices lead to different "prescrip- 
tions" . Once the other NLO corrections are included, the 
physical result should be independent of this prescription. 
However, in present day applications this is not done yet, 
and the result will depend on how the running coupling 
is implemented. The running coupling constant prescrip- 
tion that is usually used with the BK equation is the one 
derived by Balitsky [20] (our Eq. (38)), since it has been 
found to minimize a certain subset of other NLO contri- 
butions [21] and provide a slow evolution speed roughly 
consistent with HERA data. 

The few existing solutions [11, 29] of the JIMWLK 
equation with a running coupling have not used the same 
prescription. This is due to the need to decompose the 
evolution kernel into a product of two factors for the 
numerical solution of the JIMWLK equation. Instead, 
what has been used is a "square root" running coupling, 
Eq. (15). What we propose in this paper as a solution 
to overcome this technical limitation is an alternative 
form of the running coupling JIMWLK equation, given 
by Eqs. (25) and (26). 

In this paper our starting point is different from 
Refs. [19-21], where one sets out to find the correct mod- 
ification of the BK kernel, keeping the form of the BK 
equation otherwise intact. Here we want to propose an 
equation ( "rc JIMWLK" ) that keeps intact the functional 
form of the JIMWLK equation and there is no particu- 
lar reason for the resulting mean field BK approximation 
to stay exactly the same. We will, however, argue the 
BK equation derived from our proposal (Eq. (27)) is for 
practical purposes very close to the one with the Balitsky 
prescription. Its functional form is also, while different 
from rcBK, not that different from structures appear- 
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ing in the full NLO BK equation. We will also compare 
numerical simulations of both the BK equation and our 
proposed JIMWLK formulation, starting from the same 
initial condition. 

This paper is organized in the following way. We 
shall first, in Sec. II, recall the Langevin formulation 
of the JIMWLK equation existing in the literature and 
in Sec. Ill the derivation of the BK equation from it. 
We proceed in Sec. IV to rewrite the Langevin step in a 
simpler way that corresponds to the the manifestly left- 
right or "mirror" [30-32] symmetric form of JIMWLK. 
We then motivate our proposal for a running coupling 
version of the equation in Sec. V and discuss its relation 
to the rcBK equation in Sec. VI. In Sec. VII we compare 
numerical solutions of the rcBK and rcJIMWLK equa- 
tions before concluding in Sec. VIII. 



II. LANGEVIN FORMULATION OF JIMWLK 

The equation for the probability distribution of Wilson 
lines can be seen as a functional Fokker-Planck equa- 
tion [4] . This in turn can be re-expressed as a functional 
Langevin equation for the Wilson lines themselves [33], 



V x = V x (it a ) 



I ab,i j- ( \b i a 
J z 



(i) 



Here we denote two dimensional vectors by x, • ■ ■ and 
their lengths as x T . In the interest of space the coor- 
dinate arguments are written as subscripts. The Wilson 
line V is a unitary matrix in the fundamental representa- 
tion, generated by t a and i = 1, 2 is a transverse spatial 
index. The first term in Eq. (1) is proportional to a 
stochastic random noise and the second one is a deter- 
ministic "drift" term. The coefficient of the noise in the 
stochastic term is 
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which can be thought of the "square root" of the 
JIMWLK Hamiltonian. We follow here the notation of 
Ref. [11] and denote the adjoint representation Wilson 
line by U; it is related to the fundamental one by 



conditions. For this reason we will for the moment keep 
this more general notation. The noise £ = (£,it a , is 
taken to be Gaussian and local in color, transverse coor- 
dinate and rapidity y (evolution time) with (£ x (z/)i) = 
and 

(UvMyiy')*) = S ab 5 l >5$5(y y), (5) 
The deterministic term is 
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where T° is a generator of the adjoint representation, Tr 
the adjoint representation trace and we have introduced 
the kernel 



1 



(x-zr 



(7) 



The Langevin process is both derived [33] and imple- 
mented numerically [34, 35] in a discrete evolution time. 
It is therefore useful to see explicitly what happens on 
one step in rapidity from y n to y n+ i — y n + dy. The 
change in the Wilson line is 

Vjy+ dy) = V x (y) 

x exp |ii a J £x fc z l £z 1 \fd~y + &i dj/| . (8) 

Note that the rapidity delta function in Eq. (5) becomes a 
Kronecker delta divided by the magnitude of the timestep 
^(Vm ~ Vn) ~^ &m.nl which we have removed from the 
normalization of the noise to keep the power counting 
in dy more explicit. Thus the variance of the Gaussian 
noise is now 



(£xM°£y(")i) = S ab S l3 S^5„ 



(9) 



and the stochastic term in Eq. (8) is explicitly propor- 
tional to y/ dy. The step in rapidity dy is assumed to 
be small so that we can expand all the quantities on one 
timestep to the order dy, i.e. to second order in the 
stochastic term and to first order in the deterministic 
term. 



uf = 2Trt a Krv; 
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We have denoted the coordinate dependent part of the 
kernel (essentially the q —> qg splitting wavefunction in 
LC quantization) by K x _ z . For the fixed coupling con- 
tinuum theory it is simply 
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In the "square root" prescription one would modify this 
to account for the running coupling, and on a lattice one 
needs to render K x _ z consistent with periodic boundary 



III. DERIVING BK FROM JIMWLK AT FIXED 
COUPLING 

The BK equation describes the time evolution of the 
two point function or "dipolc" 



(10) 



It can be derived from JIMWLK in a mean field approx- 
imation, where four point functions of Wilson lines are 
assumed to factorize into products of two point functions; 
what matters here is, however, not this truncation of the 
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Balitsky hierarchy but the relation between the BK and 
JIMWLK kernels. Let us start by rederiving the BK 
equation directly from Eq. (8). In order to get the ra- 
pidity derivative of -D x ,y we must express D xy (y + dy) 
as -D Xj y(y) plus a correction proportional to dy. This 
is done by inserting the timestep (8) into the expression 
for D xy (y + dy) and expanding to order dy. Since we 
are, in the BK equation, interested only in the evolution 
of the mean scattering amplitude we can then take an 
expectation value over the stochastic noise £, which is 
independent of the Wilson lines from the previous steps 



in rapidity (i.e. the Langevin equation is discretized with 
the Ito interpretation). One then uses the Fierz identity 



t a Mt a = - TrM M (11) 

to remove the generators remaining from the Gaussian 
avarage (£ a £ b )i a • • • t b ■ ■ ■ = t a ■ ■ ■ t a ■ ■ ■ . To clarify the 
developments in what will follow, let us gather the terms 
in (D x , y (y + dy) - £> x , y (y))/dy here one by one: 
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^^((i*^) 2 )^ ^ ?i K ^T{^ Tr ^ ^ T ^ zVy -^ T ^ xVy 

+ ^ Tr KV X ~ Tr F x V u t/ x V y 
^TVV x V y (i e ^'^'^ = ^ /K^.K^i^Tr^Vy-^Tr^V^Tr^ 



The terms that are linear in the noise vanish in with the expectation value. Combining the equations above yields 

dy 2ir 2 J z \ \ 2 2 y ^ JV C x A c y 

' 'k 2 - S,-^ i Tr K f K i Tr V$V*VlV v + I ( K 2 _ v - S1_J\ Tr vlv. i Tr kV^V, 



2 ^ — x-z ^ x -z y^y -^'z" X ^y "'x'z'x'y i 2 ^ z -y ^z-y /7V yz 7V z y x y 

K 2 z + K 2 _ y - 2K X _ Z • K z _ y ) i, Tr V x V y 1 . (13) 



A few remarks are in place based on an inspection of 
Eq. (13). Firstly, one cannot separately modify the ker- 
nels S and K in order to discretize the theory or to in- 
troduce a running coupling. If the relation K 2 = 5" is 
not fulfilled one introduces an additional six point corre- 
lation structure into the evolution equation of the dipole. 
Although this kind of higher point functions do appear 
in the NLO BK equation, they should be corrections at 
order a 2 and not present here at the leading order, fixed 
a s level. Secondly, the stochastic term (identifiable from 
the kernel K instead of S) alone gives the correct BK 
equation virtual term, i.e. the one with only the original 



dipole on the r.h.s. The stochastic term alone, however, 
gives an incorrect kernel for the real term and the addi- 
tional sextupole term, which are then corrected by the 
addition of the deterministic term. We shall show in the 
following that the additional complication of having to 
include a deterministic term is an artefact resulting from 
writing the evolution step only in terms of a multipli- 
cation of the Wilson line from the right (right covariant 
derivative). The JIMWLK equation is in fact naturally 
left-right-symmetric, which can be made explicit by ex- 
pressing it in terms of both left and right derivatives [30- 
32]. 
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The fixed coupling BK equation, included here for ref- 
erence, can be obtained from Eq. (13) by setting S = K 2 
to get 
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K 2 +K 2 2K K 



D(x,z)D(z,y) - D(x,y) . (14) 



To obtain the BK equation for the dipole (D) one then 
closest this by the mean field approximation (DD) ps 
(D)(D). 

The "square root" running coupling prescription of 
Refs. [11, 29] was simply justified as the simplest possible 
modification of the kernel K to include a dependence of 
the coupling on the distance scale. If the only modifica- 
tion possible is that of the kernel K, one cannot impose 
a running of the coupling as a function of the "parent" 
dipole size |x — y|. In stead, one is left with the "square 
root" coupling 



d y D x>y (y) - 2 

Z7T 



" s (l x - z l) K x-z+a s (|y-z|)K y _ z 

2 v / a s (|x- z|)a s (|y - z|)K x _ z • K y _ 2 

D(x ) z)D( Z ,y)- J D(x,y)j . (15) 

The argument of the kernel K is, however, not the only 
scale in the problem that the coupling can depend on. 
Introducing an a s into the correlation function of the 
noise is another option which, as we will argue in Sec. V, 
has a nice physical interpretation. But first it is useful 
to simplify the Langevin form of the equation slightly. 

IV. SIMPLER FORM OF THE LANGEVIN STEP 

The presence of the deterministic term in Eq. (1) is 
somewhat of an inconvenience for the discussions that 
will follow. It is also a relatively expensive part of the nu- 
merical solution of the equation. As discussed in Ref. [5], 
the deterministic term is related to the need to trans- 
form "left" and "right" covariant derivatives into each 
other. In the Langevin formulation this means that it is 
a consequence of wanting to write the evolution step as a 
multiplication of the Wilson line only from one side. We 
shall now show how one can write the rapidity step with 
only a stochastic term by allowing for a multiplication of 
the Wilson line from both left and right. 

In order to derive an equation that describes the whole 
Balitsky hierarchy, not just the evolution for the dipole, 
we must start by finding out what exactly is the change 
of a single Wilson line at point x to sufficient order in 
dy. In the evolution step of a generic Wilson line opera- 



tor V, 



* y yy ■ ■ <£> Vy n one needs to know, 
in the limit dy — > 0, its change up to 0( dy)-terms; the 
higher ones vanish in the continuous rapidity limit. Now 



a 0(dy)-term can come from either the product of two 
O (\/ dy) -terms from different Wilson lines, or from a 
0(dy)-term from a single Wilson line. It follows that for 
a single Wilson line we must, in the Langevin timestep, 
keep the terms linear in the noise £, which could be con- 
tracted with the noise term coming from another Wilson 
line. One must also keep the O( dy) terms, but for them 
one can take the expectation value of the noise term, be- 
cause any correlation between them and the evolution of 
another Wilson line would be higher order in dy. Ex- 
panding Eq. (8) to order dy and taking the expectation 
value of the quadratic noise term one gets 



V x (y+ dy) = V x (y) 



i • V" s dy . 
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Here V X V Z $ Z V Z V X is simply the adjoint representation 
product t a (U x U z ) ab £ z written out in the fundamental 
representation. The terms coming from the noise and de- 
terministic contributions can be identified from the ker- 
nels K 2 and S respectively. The effect of the determin- 
istic term is to cancel U\\J X and double the other term 

t 2 

UyXJ z . Now setting S equal to K and also expressing the 
remaining adjoint Wilson lines in terms of fundamental 
ones we can write this as 

V x (y+ dy) = V x 

+ I K x _ z • (V^ z - VUzKV x ) 

7T Jz 

- [ K 2 _ Z (V Z TrKV x - N C VJ. (17) 

Up to this order in dy this is equivalent to the re- 
exponentiated form 



V x (y + dy) - cxp i I Kx _ z ■ (V m $ u V^ 



x K 



6XP {^i Kx - z ^}' (18) 



where we have again, in the order dy replaced (£ ) by 
£ , the difference being higher order in dy. 

Equation (18), to be compared to the earlier form (16), 
is the main result of this section. Now the Wilson line is 
multiplied from both the left and the right, and conse- 
quently there is no more drift term. In addition to having 
fewer separate terms to calculate, Eq. (18) is significantly 
faster to evaluate numerically, because it avoids the need 
to reconstruct the adjoint representation Wilson line in 
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Eq. (6). We stress that Eq. (18) is precisely equivalent to 
the Langevin step derived in [33] in the limit of continu- 
ous evolution time dy 0. The difference comes only in 

terms at order O (^/dy 3 ^j , which are neglected both here 

and in the derivation of the original equation in Ref . [33] . 

Note that the multiplication by the Wilson line at the 
coordinate z from the left is a choice related to the form of 
the original choice to use a right derivative in (1). One 
could easily move it to the right by defining a rotated 
noise 

L=U a = V z £ z Vl (19) 

It is easy to see that the distribution of £ z is the same as 
that of £ z , namely 

= 6 ab ^6?l y . (20) 
In terms of £ z the timestep (18) is 

V x (y + dy) = exp {-i^^ jf Kx - Z ' *>} 

x y x exp|z^^^K x _ z .(^ z K)}. (21) 



V. SETTING THE SCALE OF THE COUPLING 

Recall that the change in a general correlator of n Wil- 
son lines in one infinitesimal step in rapidity is obtained 
by changing all the Wilson lines by Eq. (18), develop- 
ing to order £ (i.e. to order dy) and taking the ex- 
pectation values with the probability distribution of the 
noise £. Physically keeping only the quadratic order in 
£ corresponds to the fact that JIMWLK is a LO evo- 
lution equation, derived by considering the emission of 
only one gluon. The interpretation of the contractions 
between two £'s is that one corresponds to the emission 
of the gluon in the amplitude and the other one in the 
complex conjugate. In the formulation that includes the 
deterministic term this interpretation is not as straight- 
forward, since in Eq. (1) one is summing a noise term 
corresponding to a gluon emission amplitude and the de- 
terministic term which is an emission probability, i.e. am- 
plitude squared. 

In order to have a better picture of the physics con- 
tained in one evolution step we shall now describe, adapt- 
ing the "simple derivation of the JIMWLK equation" 
of Ref. [5], the contractions of the £'s in Eq. (18) in a 
diagrammatic notation. We want to apply the evolution 
step Eq. (18) to calculate the change of a generic Wilson 
line operator V% ® • • • <8> V-^ n <S> V yi ® • • • ® V y in an in- 
finitesimal step in rapidity. We now have to distinguish 
between several different cases of how these contractions 
are done, drawn in Figs. 1 and 2: 



• The square of the noise multiplying the Wilson line 
from the left corresponds to the emission and rea- 
sorption of a gluon before interaction with the tar- 
get. The noise term £ z is color rotated into V Z £ Z V Z 
in Eq. (18), but in the square the color rotation can- 
cels. Thus we can express this contribution as the 
left diagram in Fig. 1. 

• Similarly we can take the square of the £ multiply- 
ing the Wilson line from the right. Now the noise is 
not color rotated, i.e. this corresponds to emission 
from and absorption by the Wilson line after the 
target, shown in the middle diagram of Fig. 1. 

• We can also contract the £ multiplying the Wilson 
line from the left with the one multiplying it from 
the right. Physically this is equivalent to the in- 
terference between gluon emission before and after 
the interaction with the target. In the diagram- 
matic notation this is expressed as a gluon being 
emitted before the target, propagating though it, 
and being reabsorbed, drawn on the right in Fig. 1. 

• One can also have a contraction between the £'s 
from the step of one Wilson line at coordinate x 
with another one at y. If both £'s multiply the 
Wilson lines from the left (or the hermitean con- 
jugates from the right), they correspond to the di- 
agram where the gluon is emitted and absorbed 
before the target, the color rotations again cancel. 
This contribution is shown on the left of Fig. 2. 

• If the contracted £'s are the ones to the right of 
the Wilson lines, we have an interference between 
gluon emission from two different Wilson lines after 
the target, as shown in the middle of Fig. 2. 

• Finally we can have an interference between the 
emission of a gluon from the Wilson line at x before 
the collision and from another one at y after the 
collision, as shown on the right of Fig. 2. This 
corresponds to the contraction between one £ on 
the left and another one on the right of a Wilson 
line. 

If now £ u corresponds to the emission of a gluon at 
coordinate u and the other £ v that it is to be con- 
tracted with is interpreted as the absorption of the same 
gluon, what is the physical meaning of the delta func- 
tion (5 2 (u — v) in the correlation function of the noise? 
The answer to this question that we propose in this pa- 
per lies in realizing that in fact the "absorption" of the 
gluon actually corresponds to the complex conjugate of 
the amplitude for gluon emission. Thus if the transverse 
momentum of the emitted gluon is k, the amplitude is 
~ e lk u and the complex conjugate ~ e * k v . The delta 
function S 2 (u — v) is then a result of integrating over all 
possible transverse momenta of the emitted gluon in one 



FIG. 1: Virtual diagrams. 
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FIG. 2: Real diagrams. 



step of the evolution: 



<T(u- 



now 



d 2 k 

(27T) 2 



ik-(u— v) 



(22) 



The interpretation of the diagrams of Figs. 1 and 2 in 
terms of the transverse momentum of the gluon is shown 
in Figs. 3 and 4. Note that all the momenta k are taken 
after the interaction with the target, not before. This 
corresponds to the choice made to use the noise £ and 
not the rotated one £. The integration over the final 
state gluon momentum has to be done because we wish to 
include all the degrees of freedom in the rapidity interval 
from y to y + dy into an effective theory description of 
the target at the new rapidity scale dy. 

What then happens if we wish to include the running 
coupling constant? We will here rely on the general ar- 
gument in gauge theory that the beta function can be 
computed by considering only corrections to the gluon 
propagator. This is in fact the property that is used 
in [19, 20] to derive the running coupling constant in BK 
which is done by summing one loop corrections to the 
gluon propagator. If the running coupling corrections 
can be derived from just the gluon propagator, it is nat- 
ural that the scale of the running coupling should be the 
momentum of the emitted gluon. This leads directly to 
our proposal for implementing the running coupling in 
the JIWMLK equation, which is to regard the coupling 
constant as a property of the correlator of the noise. We 
then propose to simply replace the fixed coupling corre- 
lator 



S ab aJ l] 
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(27T) 2 



ik-(x-y) 



(23) 



V X (V + dy) = cxp j-i^ j K x _ u • (V uVu V^ 

.Vdy- 



x V x exp < i 



K x v • »7 V 



(25) 



with r\ — (rjit a , i]%t a ). Note that since the noise correla- 
tor is not local in transverse coordinate, we have broken 
the left-right symmetry of the fixed coupling equation. 
The physical interpretation of this is that the symmetry 
corresponds to time reversal, and using the momentum of 
a gluon in the final state after the target (as opposed to 
an initial state one) breaks the time reversal symmetry. 
While this loss of symmetry in inacsthctic, it seems in- 
evitable at some point at least if one wants to generalize 
the evolution to more exclusive observables (as discussed 
e.g. in Ref. [36]). 

Due to the FFT algorithm [37] the Fourier-transform 
is a relatively straightforward and inexpensive procedure 
on the lattice. Thus implementing the correlator (24) 
induces only a minor modification to the numerical algo- 
rithm [34] used to solve the JIMWLK equation. In the 
following we shall denote the new noise correlator by 



d 2 k 

(27T) 2 



ik-(x-y) 



a s (k). 



(26) 



We emphasize that the notation 5 x _ y does not mean 
the coupling evaluated at the momentum scale l/|x — y|, 
but is indeed the Fourier transform of the coupling at 
the scale k T . Since the coupling is a very smooth, log- 
arithmically varying function of fop, the transform a x _ y 
is very sharply peaked around x = y. 



by 



(r,rVy' j ) = S ab 6 ij j -|^ e ik -(*-r)a 8 (k), (24) 



(27T) 



in terms of which the rapidity step of the Wilson line is 



VI. RECOVERING THE BK EQUATION 

Let us now calculate the change in the dipole operator 
during one rapidity step of the evolution. This is done by 
inserting Eq. (25) into the expression for D xy (y + dy), 
developing in powers of the noise term up to 0(r) 2 ) and 
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FIG. 3: Virtual diagrams in Fig. 1 rewritten with the explicit momentum argument used as the scale of the running coupling. 
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FIG. 4: Real diagrams in Fig. 2 rewritten with the explicit momentum argument used as the scale of the running coupling. 



taking the expectation values with the correlator (26). The result is 



V f „ f 

- 2 / — (K X _ U ' K X _ V + K y _ u ' K y _ v — 2K X _ U • K y _ v ) 

7T J u ,v [ 
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X - 
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r 



(27) 



where the quadrupole operator is 



1 



(28) 



From this form it is easy to check that in the fixed cou- 
pling case 5 U _ V = a s S 2 (u — v) we recover the usual BK 
equation. 

Now the evolution equation for the two point func- 
tion does not have the same functional form as in the 
fixed coupling case. Our definition of the scale of the 
running coupling as the momentum of the emitted final 
state gluon makes the position of the gluon in transverse 



coordinate space different in the amplitude and the com- 
plex conjugate, resulting in a more complicated the struc- 
ture for the equation. From the point of view of the full 
set of NLO corrections to the BK equation [27] this is 
not unnatural. The NLO BK equation also necessarily 
involves higher point dipole operators, and the defini- 
tion of what subset of these corrections one chooses to 
call "running coupling" corrections is not unique. What 
is done in the derivations of the running coupling BK 
equation [19, 20] is that one sets out to calculate a cer- 
tain subset of these NLO corrections, namely those that 
only modify the BK kernel, but leave the Wilson line 



operator part of the BK equation intact. This partic- 
ular set of NLO corrections is not necessarily one that 
would leave the functional form of the JIMWLK equa- 
tion intact. What we are proposing here, on the other 
hand, is running coupling JIMWLK equation, namely an 
evolution equation which includes running coupling cor- 
rections, but leaves the functional form of the JIMWLK 
equation intact. There is no particular reason why this 
equation should be equivalent to a particular form of the 
running coupling BK equation, and indeed it is not. How- 
ever, as we will show in the following, it does in fact lead 
to the scale of the coupling being determined by para- 
metrically the same length scale in the parts of phase 
space that dominate the evolution. 

In order to see the relation to rcBK equations in the 
literature in more detail we will now try to connect our 
Eq. (27), where the coupling is a function of a momentum 
scale to one where it depends on the inverse of a trans- 
verse coordinate distance. We shall here use the same 
notation for both, a s (k) and a s (r), while deferring the 
explicit expressions to Sec. VII. We will also not attempt 
to perform the Fourier-transform from momentum to co- 
ordinate space precisely, but merely attempt to identify 
the dominant momentum scales in parametrically differ- 
ent parts of the phase space. 

Let us denote z = (u + v)/2 and Az = (u — v). Using 
the fact that 5 U _ V is very peaked around u = v we then 
approximate u w v ss z (i.e. set Az = 0) in the part of 
(27) containing the Wilson line operators 



-^x,u-^u,y ~1" -^x,v-^v,y -^x,y ^v,uQx,v,u,y 

AcAy-Ac.y, (29) 



U— V — z 



which is the same form as in the BK equation. 

Turning then to the kernel, we can use the momentum 
representation of K 



K x = i(2n) 



d 2 k e"' kx 
(2vr) 2 k 2 



(30) 



to write 



/ a Az K x _ u • K y _ v 

JAz 



(27T) 



p,q 
1 



P + Q ^ P • Q -ip x+iq y+l(p-q) z 
2 2 e 

p q 



„ , a s (P)e 
(2irfJ P M 



iP-(x-y)-Mk-(z-5±£) 

(P + k/2)-(P-k/2) 
(P + k/2) 2 (P-k/2) 2 



(31) 



where we have changed variables from p, q to P = (p + 
q)/2 and k = p — q. Using this form we can write the 



kernel of Eq. (27) as 



«Az (K x _ u ■ K x 



Az 



y— u ' ^ v v 



2K X _ U • Ky_ v ) 



1 



(P + k/2)-(P-k/2) 
(2^) 2 7 P , k ; (P + k/2) 2 (P-k/2) 2 

e ik-(z-x) e ik-(z-y) _ 2 e -«P-(x-y)+ik-(z-^jp 



(32) 



We now look separately at two different kinematical 
regimes. 

a. Small daughter dipole Let us first assume that 
one of the daughter dipoles is smaller, e.g. |x — z| <C 
|y — z| w |x — y|. Now the dominant term in the P and 
k-integral is the one that oscillates most slowly, namely 
as e ,k '( z_x ). Keeping only this term we can return to the 
variables p and q and write it as 



1 



(27T) 



P • q -t(p-q).(x- g ) 



2 2 ' 

p q 



(33) 



The exponential factor now makes both p T and of the 
order l/|x — z|. We thus approximate a s ((p + q)/2) « 
a s (|x— z|), leading us to the estimate 



a s (x-z) 



|x-z| < y-z 



(x-zr 

for the kernel in this limit, and similarly 

« s (y - z 



z < x - 



(34) 



(35) 



for the other small daughter dipole limit. This is almost 
an inevitable result, which can also be seen noting that 
already at fixed coupling the kernel in this kinematical 
limit is dominated by the single term K 2 z or K 2 _ z , 
and the coupling must depend on the only distance scale 
available for it, namely the daughter dipole size. 

b. Small parent dipole Assume then that the parent 
dipole is smaller, e.g. |x — y| <C |x — z| ss |y — z|. We 
can then approximate 

g ik-(z-x) _j_ e ik-(z-y) _ 2 e -<P-(x-y)+»k-(»- s y x ) 

s , 2e jk ( z -^ ) (l- e - iP(x - y )), (36) 

from which we immediately see that the scale of the cou- 
pling P is given by the parent dipole size x — y. The 
physical intepretation of this is that the emission of the 
additional gluon happens coherently from the dipole (as 
evidenced by the fact that the kernel vanishes at x = y, 
reflecting color neutrality), thus the parent dipole deter- 
mines the scale of the coupling. After changing the scale 
of the running coupling to x — y we can then pull it out- 
side of the P, k-integral (32) and restore the exponential 
functions to the original form on the l.h.s. of Eq. (36). 
Doing this we get the estimate 



a 8 (x-y) 



(x 



(x- 



x-y < x-z 



(37) 
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for the kernel. 

It is instructive to compare the estimate (34), (35) and 
(37) to the "Balitsky" rcBK equation, which reads, in our 
notation: 



N 

d y D x . y (y) = 



g s (x-z)a s (x-y) 2 

" s (y- z ) 



« s (y- z)« s (x-y) K 2 
a s {z-z) y ~ z 



2a s (x-y)K x _ z -K y _ 2 

x £>( X) z)D( Z ,y)-D(x,y) . (38) 



One immediately recognizes that the scale of the coupling 
in the different kincmatical limits is the same. 

As a summary, showing the approximate equivalence of 
our rcJIMWLK equation, defined by Eqs. (25) and (26), 
to the "Balitsky" rcBK equation (38) commonly used in 
phenomenology involved the following steps. We first de- 
rived an equation for the two point function Eq. (27) that 
involved 4 transverse coordinates and a six-point function 
DQ. We subsequently took the limit of only 3 coordi- 
nates, in stead of 4, for the Wilson line operators, reduc- 
ing them to the BK form. We then identified the domi- 
nant momentum scales for the coupling in the kernel of 
Eq. (27) in the limits of small and large daughter dipoles 
compared to the parent dipole, and interchanged the mo- 
mentum and the corresponding coordinate distance as 
the argument of the running coupling. In addition, our 
Eqs. (27) and (38) have been written as operator equa- 
tions corresponding to the first equation of the infinite 
Balitsky hierarchy. While the JIMWLK equation solves 
the whole hierarchy simultaneously, a solution of just the 
BK equation requires one to truncate the hierarchy by re- 
placing the operators D by their expectation values. To 
quantify the effect of these approximations we shall now 
proceed to a comparison between numerical evaluations 
of the BK and JIMWLK equations. 



VII. NUMERICAL COMPARISON 




QCD 



FIG. 5: Evolution of the scattering amplitude with the 
JIMWLK equation where the scale of the running coupling 
is set in the noise term, Eq. (24), compared to BK evolution 
with the "Balitsky" running coupling prescription, Eq. (38). 
The evolution is started from the same initial condition and 
the lines show the amplitude at intervals of 2 units in rapidity. 




QCD 



FIG. 6: Evolution of the scattering amplitude N = 1 — (D T ) 
with the rcJIMWLK equation proposed in this paper (same 
data as in Fig. 5), compared to the "square root" prescription 
(15) used in [11, 29]. The evolution is started from the same 
initial condition and the lines show the amplitude at intervals 
of 4 units in rapidity. 



The numerical algorithm for solving the JIMWLK 
equation is discussed in detail in [34], and the modifica- 
tions required to adapt it to our version of the JIMWLK 
equation, Eqs. (25) and (26), are straightforward. We 
construct the initial condition Wilson line configurations 
as described in Ref. [38], using a 1024 2 transverse lattice, 
with the MV model parameter g /j,L = 31 and the longi- 
tudinal extent in the MV model discretized in N y = 100 
steps, leading to an initial fundamental representation 
saturation scale Q S L ss 22, where L is the size of the 
lattice in physical units. The infrared behavior of the 



coupling is regularized as 



a s (k) = 




with (3 = 11 - 2iV F /3, where we take N F = 3. The 
parameter c = 0.2 regulates the sharpness of the freezing 
of the coupling, which happens at /i L = 15. The value 
of the QCD scale is taken as Aq CD L = 6, leading to a 
value a = 0.7619 for the coupling in the infrared. In the 
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FIG. 7: The momentum space correlator, (41), from the 
rcJIMWLK equation proposed in this paper compared to the 
"square root" prescription (with the same parameters as used 
in Fig. 6). The correlator is shown for the initial condition, 
for which we use the MV model, and for y = 10 and y — 7.8 
units of evolution for the two running coupling versions; these 
correspond to approximately the same value of Q s . Note that 
the end of the curves at large k T corresponds to the lattice 
ultraviolet cutoffs, and the shape is sensitive to lattice effects 
well below that. 
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FIG. 8: Evolution speed of the saturation scale as a function 
of Q s /Aqcd- Shown are the noise term (solid line), square 
root (dotted line) running coupling JIMWLK equations and 
the Balitsky prescription BK equation (dashed line), starting 
from the MV model at two different initial values of Q s . The 
thin lines have the freezing of the coupling implmented more 
smoothly, with the value c = 1.5 instead of c = 0.2 as used 
elsewhere in this paper. 



more extensive study in Ref. [29] the evolution speed was 
found to be very insensitive to the lattice parameters and 
to chq while naturally dependent on the physical ratio of 
scales Q s /Aq CD . 

When the coupling is needed as a function of the trans- 
verse coordinate, as in the BK equation or the ^/cig- 
prescription for JIMWKL that we study as a comparison, 



we use 

«s( r ) 



4tt 



/31n- 



»QCD 



4e 



r A 



QCD 



(40) 



Note the constant factor 4e 

-2 7E 



-2 7e 



1.26 in the identifi- 



cation fcy ~ 4e~ Z7E /r T , which is taken from the explicit 
Fourier-transform of the kernel in Refs. [19, 39]. In Ap- 
pendix A we show a comparison between the coordinate 
and momentum versions of the running coupling, verify- 
ing that this is indeed the correct correspondence. Phe- 
nomenological studies typically [22, 23] find that the co- 
ordinate space scale needs to be taken as larger, ~ 20/r^ 
for Aq CD > 200MeV, unless some other corrections are 
applied that slow down the evolution, such as in Ref. [25]. 
Thus, in spite of the slowness of the rcBK evolution with 
the Balitsky prescription compared to fixed coupling, the 
data seems to prefer a slower evolution still than it would 
predict. 

Figure 5 compares the evolution of the dipole using the 
rcJIMWLK equation (Eqs. (25) and (26) with (39) for the 
coupling) with the Balitsky prescription rcBK equation 
(38) with the coupling (40) over 20 units in rapidity. We 
see that the Balitsky prescription leads still to a slower 
evolution than our proposal, in spite of the similarity 
in parametrically different parts of the phase space dis- 
cussed above. Discretization effects that also affect the 
comparison are discussed in Appendix A. The shape of 
the correlator as a function of r T is steeper in BK than 
in JIMWLK, for which we do not see a compelling ex- 
planation. 

Figure 6 compares the rcJIMWLK evolution using our 
proposal to the "square root" coupling (15) used in pre- 
vious rcJIMWLK calculations [11, 29]. One can see that 
the the evolution speed is slowed down significantly from 
the " square root" prescription" . While the shape of the 
front seems somewhat different there is no clear sign of a 
different power law in the momentum space correlator 



C(k) 



,2 ik-j 

d re 



D(r T ) 



(41) 



This is shown in Fig. 7 as a function of the scaling vari- 
able k T /Q s for both coupling constant prescriptions at a 
similar value of Q s (which, due to the different evolution 
speed, correspond to slightly different values of rapidity). 
More detailed statements on the geometric scaling behav- 
ior at large k T would require significantly larger lattices 
to push out the ~ 1/a ultraviolet cutoff further. 

Figure 8 further illustrates the difference in evolution 
speeds between the prescriptions. The evolution speed is 
defined as 



A 



d\nQj 
dy 



(42) 



with the saturation scale Q s solved from the condition 



D 



2 

r T 



2/q; 



-1/2 



(43) 
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In addition to the initial condition used in Fig. 6, evolved 
for 20 units in y, Fig. 8 also shows the result of a calcula- 
tion started with an initial Q s twice as large. The evolu- 
tion speed shows the typical behavior also seen with the 
BK equation (see e.g. [25, 40]), with a slower evolution in 
the initial transient region when the dipole cross section 
changes from the initial condition to the pseudoscaling 
form, followed by a universal 1/^/y-dependence charac- 
teristic of rcBK evolution. For comparison Fig. 8 also 
shows the evolution speed calculated with a smoother 
freezing of the coupling in the infrared, implemented by 
choosing c = 1.5 in Eqs. (39) and (40). The depen- 
dence on the freezing prescription, which is not present 
for larger initial Q s /Aq CD (not shown on the figure) is, 
interestingly, different for the BK and JIMWLK calcula- 
tions. This initial behavior of course also depends on the 
functional form of the initial condition. At large rapidi- 
ties, the evolution speed in JIMWLK drops below that 
in BK due to the fact that Q s grows close to the lattice 
ultraviolet cutoff; see the discussion in Appendix A. 

We remind the reader that the typical phenomenolog- 
ically favored initial saturation scale for the evolution is 
Q s ~ 1 GeV and the evolution speed resulting from fits 
to HERA data (e.g. [22, 23, 25, 41-43] ) is A - 0.2 . . . 0.3. 
The "square root" coupling could be made to fit experi- 
mental data by adjusting the coefficient under the loga- 
rithm in the running coupling as ln[4C 2 /(r 2 Aq CD )] effec- 
tively taking Aq CD as a fit parameter. Doing so would, 
however, lead to an unnaturally large value for C, i.e. 
small value for the effective Aq CD . With our rc JIMWLK 
prescription, as with the Balitsky rcBK, this problem is 
much less severe, although still present to some degree. 



VIII. CONCLUSIONS 

Until now running coupling solutions of the JIMWLK 
equation have been limited by the belief that since the 
Langevin formulation is written in terms of the square 
root of the JIMWLK kernel, one is limited to consider- 
ing only prescriptions that modify this square root. In 
particular, since the square root only depends on the size 
of the daughter dipole, it has been impossible to include 
a dependence on the "parent" dipole size in the coupling. 
We have in this paper pointed out that there is in fact 
another scale also in the Langevin formulation that the 
running coupling can depend on, namely the scale in the 
random noise term. We have proposed to use this scale as 
the argument of the running coupling, with the Langevin 
timestep of the JIMWLK equation given by Eqs. (25) and 
(26). Choosing this as the scale of the coupling has sev- 
eral advantages. Firstly it corresponds naturally to the 
general argument that the scale of the coupling should be 
the momentum of the emitted gluon. Secondly we have 
shown analytically how it leads to the typical scale of 
the coupling being taken as the smallest one of the three 
relevant dipole sizes (the "parent" and two "daughter" 
dipoles) present in a dipole splitting. This is very similar 
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QCD 

FIG. 9: Evolution of the scattering amplitude N = 1 — 
(D T ) with the "square root" coupling rcJIMWLK equation, 
with coupling evaluated at the scale in the momentum 
space kernel, and at the scale 4e _27E /ry in the position space 
kernel. 



to the "Balitsky" prescription for the running coupling 
in the BK equation. We have shown numerically that 
the evolution speed of the saturation scale in our pro- 
posed running coupling JIMLWK equation is slower than 
the "square root" coupling, although not fully consistent 
with the Balitsky prescription for BK, or with fits of the 
data [22, 23]. 
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Appendix A: Lattice effects 

As test of the equivalence of the scales in the mo- 
mentum space coupling and the scale 4e 27e jr% in coor- 
dinate space we show in Fig. 9 a comparison between a 
solution of the rcJIMWLK equation using the two cou- 
plings. This can be done strainghtforwardly when one 
uses the "square root" coupling, where one can easily 
multply either the kernel K x or its Fourier-transform 
with y/a^. Thus the curves in Fig. 9 correspond to evo- 
lution using the kernels 

AWK r , (Al) 
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FIG. 10: Evolution of the scattering amplitude N = 1 — (D r ) 
with the "square root" rcJIMWLK equation (dotted line) and 
with the BK equation (thinner dashed line) using the same 
"square root" running coupling. The initial condition is the 
same as in Fig. 6 and the amplitude is shown every 4 units in 
y. Also shown are JIMWLK simulations on a smaller 512 2 - 
lattice with the physical lattice size reduced ( "small L" , thick 
dashed line) or lattice spacing increased ( "large a" , thin solid 
line) by a factor 2. 
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FIG. 11: Evolution speed of the saturation scale A with 
the "square root" coupling. Shown are the results from the 
BK and JIMWLK simulations with the same kernel, using 
infrared regularizations of the running coupling with diferent 
values of the parameter c. 



labeled r x faZ, and 



z(2tt) 5 



d 2 k e _ik ' : 

~ W a s (k) jj- 

(2tt) 2 k 2 



(A2) 



labeled k^/a^, with a s (x) and a s (k) given by Eqs. (39) 
and (40). They are seen to agree remarkably well, jus- 



tifying our use of the constant 4e 27e in comparing the 
momentum and position space formulations. 

The BK equation is solved on a logarithmic grid in r T , 
with the values of the scattering amplitude between the 
grid points obtained from interpolation. This allows one 
in practice to study as small or large r T -values as desired. 
In the case of JIMWLK, on the other hand, one works on 
a fixed linear lattice with cutoffs ~ 1/L in the infrared 
and \ ja in the ultraviolet. One expects large lattice 
effects when Q s ~ 1/L and when Q s ~ 1/et, and since Q s 
changes exponentially during the evolution, this restricts 
the accessible rapidity range. 

To illustrate the numerical uncertainties in our calcu- 
lation we compare in Fig. 10 the rcJIMWLK calculation 
with a square root running coupling, shown also in Fig. 6, 
to the BK equation with the same running coupling pre- 
scription. One can see that the solutions agree very well 
for smaller rapidities, but start to deviate more later in 
the evolution. Also shown is a rcJIMWLK calculation 
with the same initial 5 2 A t /Aq CD (approximately the same 

Qs/Aqcd) performed on a smaller 512 2 -lattice in such a 
way that either the IR ~ 1/L or UV ~ 1/a cutoff is dif- 
ferent. We see that for the initial condition Q s is so small 
that changing the IR cutoff starts to have an effect. More 
importantly, late in the evolution the solution becomes 
sensitive to the UV cutoff; an extrapolation to a — » 
would bring the JIMWLK result closer to BK. 

A further comparison of the evolution speed in the 
BK and rcJIMWLK codes is shown in Fig. 11. It can 
be seen that in BK the evolution is slightly faster, and 
increadingly so at large Q s , when the lattice ultravio- 
let cutoff in the JIMWLK code slows down the evolu- 
tion. In a JIMWLK simulation closer to the continuum 
limit the difference would be smaller. The figure also 
demostrates the dependence on the parameter c control- 
ling the smoothness of the infrared regularization of the 
coupling. As was seen previously in Fig. 8, the evo- 
lution speed is sensitive to the regularization for small 

In addition, the BK and JIMWLK equations are not 
equivalent even in the continuum, because of the mean 
field approximation. Making a precise statement about 
this difference would require a more systematical extrap- 
olation of the lattice calculation to the continuum and in- 
finite volume limits than is done here. For fixed coupling 
the difference in the evolution speed between the BK and 
JIMWLK equations has been studied in Ref. [35], where 
JIMWLK evolution was found to be slower by few per 
cent, while an even smaller difference was seen for run- 
ning coupling in Ref. [44]. 
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